Home

workflowr

• Summary
• Checks
• Past versions

Last updated: 2024-07-17

Checks: 6 1

Knit directory: CocaineIVSA_RNASequencing_mDANeurons/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

---

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20240716)

The command set.seed(20240716) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: absolute

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute	relative
/home/avm27/Documents/Raw_Sequencing_Data/CIVSA2/CocaineIVSA_RNASequencing_mDANeurons	.

Repository version: 292e3b5

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 292e3b5. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:

Ignored files:
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  VennDiagram.2024-07-17_07-40-38.log
    Untracked:  VennDiagram.2024-07-17_07-40-39.log
    Untracked:  VennDiagram.2024-07-17_07-58-38.log
    Untracked:  VennDiagram.2024-07-17_07-58-44.log
    Untracked:  VennDiagram.2024-07-17_07-59-00.log
    Untracked:  VennDiagram.2024-07-17_08-04-55.log
    Untracked:  VennDiagram.2024-07-17_08-11-08.log
    Untracked:  VennDiagram.2024-07-17_08-22-35.log
    Untracked:  VennDiagram.2024-07-17_08-36-25.log
    Untracked:  VennDiagram.2024-07-17_08-47-16.log
    Untracked:  data/CIVSA_normalized_counts_DESeq2.txt
    Untracked:  data/gene_count_matrix_CIVSA.csv
    Untracked:  data/sample_info_CIVSA.csv
    Untracked:  output/GO_all_Combined_Comparisons.txt
    Untracked:  output/GO_all_Combined_Comparisons_Reduced.txt
    Untracked:  output/Gene_Matrix_EntrezIDs_BACKGROUND.txt
    Untracked:  output/Gene_Matrix_EntrezIDs_CRVvsFT.txt
    Untracked:  output/Gene_Matrix_EntrezIDs_CRVvsMN.txt
    Untracked:  output/Gene_Matrix_EntrezIDs_MNvsFT.txt
    Untracked:  output/Sig_Res_Craving_vs_FoodTrained.csv
    Untracked:  output/Sig_Res_Craving_vs_Maintenance.csv
    Untracked:  output/Sig_Res_Maintenance_vs_FoodTrained.csv
    Untracked:  plots/
    Untracked:  results/

Unstaged changes:
    Modified:   analysis/_site.yml

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

---

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/index.Rmd) and HTML (docs/index.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	292e3b5	avm27	2024-07-17	Testing
html	12fbe39	avm27	2024-07-17	Build site.
Rmd	479c9dd	avm27	2024-07-17	Testing
html	25a76da	avm27	2024-07-17	Build site.
Rmd	89ef098	avm27	2024-07-17	Testing
html	cb49269	avm27	2024-07-17	Build site.
Rmd	ca44efa	avm27	2024-07-17	Testing
html	21b70e3	avm27	2024-07-17	Build site.
Rmd	ddb0db7	avm27	2024-07-17	Testing
html	a1c090c	avm27	2024-07-17	Build site.
Rmd	adbd778	avm27	2024-07-17	Testing
html	dc04b15	avm27	2024-07-17	Build site.
Rmd	65a97a6	avm27	2024-07-17	Testing
html	e58ce5e	avm27	2024-07-16	Build site.
Rmd	ad8ae96	avm27	2024-07-16	Start workflowr project.

---

Cocaine craving drives repressive epigenetic differential gene expression in ventral tegmental area dopamine neurons

Tate A. Pollock1,2, Alexander V. Margetts1-3, Samara J. Vilca1-2 & Luis M. Tuesta1-3*

1 Department of Psychiatry & Behavioral Sciences 2 Center for Therapeutic Innovation 3 Sylvester Comprehensive Cancer Center University of Miami Miller School of Medicine, Miami, FL 33136 * Corresponding author (ltuesta@miami.edu)

Abstract

Dopamine (DA) signaling plays an essential role in reward valence attribution and in encoding the reinforcing properties of natural and artificial rewards. The adaptive responses from midbrain dopamine neurons to artificial rewards such as drugs of abuse are therefore important for understanding the development of substance use disorders. Drug-induced changes in gene expression are one such adaptation that can determine the activity of dopamine signaling in projection regions of the brain reward system. One of the major challenges to obtaining this understanding is the inherent cellular heterogeneity in the brain, where each neuron population can be defined by a distinct transcriptional profile. To bridge this gap, we have adapted a virus-based method for labeling and capture of dopamine nuclei, coupled with nuclear RNA-sequencing and bioinformatics analyses, to study the transcriptional adaptations, specifically, of dopamine neurons in the ventral tegmental area (VTA) during cocaine taking and cocaine craving using a mouse model of cocaine intravenous self-administration (IVSA). Our results show significant changes in gene expression across non-drug operant training, cocaine taking, and cocaine craving, highlighted by an enrichment in expression of repressive epigenetic modifying enzymes during cocaine craving. Immunohistochemical validation further revealed an increase of H3K9me3 deposition in DA neurons during cocaine craving. These results demonstrate that cocaine-induced transcriptional adaptations in dopamine neurons vary by phase of self-administration and suggest that such an approach may be useful in identifying relevant phase-specific molecular targets to alter the behavioral course of substance use disorders.

Differential Gene Expression Analysis with DESeq2

0.1 Setup

Building DESeq2 input files from gene count matrices output by StringTie following removal of samples that are outliers.

Running DeSeq2 analysis with factor levels of FoodTrained vs Maintenance, Craving vs Maintenance and Craving vs Foodtrained

1 Plot Estimates and Gather Results

1.1 Gather Results

These results include overview for each result file and dispersion estimates based on count values.

1.1.1 Craving vs FoodTraining

Craving vs FoodTrained CIVSA Results Overview

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
Mrpl15|Mrpl15	2162.26767	0.3962276	0.3468265	1.1424376	0.2532722	0.9784152
Gm39587|Gm39587	20.30129	-2.1136742	3.3911278	-0.6232954	0.5330904	0.9997970
Lypla1|Lypla1	994.90242	0.0198799	0.8799100	0.0225931	0.9819748	0.9997970
Xkr4|Xkr4	4459.42848	-0.0337912	0.2715766	-0.1244260	0.9009780	0.9997970
Gm39585|Gm39585	5797.38480	0.0861245	0.2414588	0.3566843	0.7213282	0.9997970
Tcea1|Tcea1	2657.79689	-0.4037632	0.3330614	-1.2122786	0.2254058	0.9485680
Rgs20|Rgs20	422.79650	-0.1889366	0.7619761	-0.2479561	0.8041684	0.9997970
Gm16041|Gm16041	21.30106	7.5945314	3.6799019	2.0637864	0.0390380	0.6410945
Atp6v1h|Atp6v1h	10659.10495	-0.4956230	0.2967682	-1.6700675	0.0949060	0.8425404
LOC108167788|LOC108167788	11.33585	-0.8220772	3.4983015	-0.2349933	0.8142140	0.9997970
Npbwr1|Npbwr1	924.87540	-0.3079659	0.6636912	-0.4640199	0.6426335	0.9997970
Oprk1|Oprk1	15661.74234	-0.1442457	0.3231485	-0.4463757	0.6553259	0.9997970
St18|St18	279.69289	-2.0082890	0.9250037	-2.1711146	0.0299225	0.5845690
Rb1cc1|Rb1cc1	12217.86115	-0.0796796	0.2844821	-0.2800866	0.7794110	0.9997970
4732440D04Rik|4732440D04Rik	5369.16278	-0.4918180	0.2427412	-2.0261004	0.0427545	0.6632640
Gm26901|Gm26901	163.46252	-0.2293404	0.9562622	-0.2398300	0.8104620	0.9997970
Pcmtd1|Pcmtd1	4321.05905	0.2034726	0.2645104	0.7692425	0.4417494	0.9997970
Rrs1|Rrs1	625.86245	-0.4850563	0.6605362	-0.7343371	0.4627433	0.9997970
Mybl1|Mybl1	184.22714	1.4451732	1.2066675	1.1976566	0.2310507	0.9561403
Adhfe1|Adhfe1	452.68051	-0.1536903	0.7879557	-0.1950494	0.8453543	0.9997970

out of 23922 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 250, 1% LFC < 0 (down) : 48, 0.2% outliers [1] : 2272, 9.5% low counts [2] : 1352, 5.7% (mean count < 2) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

1.1.2 Maintenance vs FoodTraining

Craving vs FoodTrained CIVSA Results Overview

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
Mrpl15|Mrpl15	2162.26767	0.6589546	0.3707556	1.7773286	0.0755142	0.3782100
Gm39587|Gm39587	20.30129	-1.7619611	3.6258386	-0.4859458	0.6270056	0.9345721
Lypla1|Lypla1	994.90242	0.0914200	0.9407122	0.0971817	0.9225821	0.9949371
Xkr4|Xkr4	4459.42848	0.6341912	0.2902594	2.1849118	0.0288953	0.2180911
Gm39585|Gm39585	5797.38480	1.1798031	0.2580070	4.5727554	0.0000048	0.0001973
Tcea1|Tcea1	2657.79689	0.1319203	0.3560044	0.3705581	0.7109667	0.9583981
Rgs20|Rgs20	422.79650	0.3482376	0.8144256	0.4275868	0.6689520	0.9469053
Gm16041|Gm16041	21.30106	7.1396504	3.9105846	1.8257246	0.0678918	0.3580246
Atp6v1h|Atp6v1h	10659.10495	0.1185235	0.3172387	0.3736100	0.7086945	0.9576220
LOC108167788|LOC108167788	11.33585	-0.5355535	3.7403736	-0.1431818	0.8861466	0.9949371
Npbwr1|Npbwr1	924.87540	-0.6449281	0.7097813	-0.9086292	0.3635459	0.7811879
Oprk1|Oprk1	15661.74234	0.2927281	0.3454541	0.8473720	0.3967878	0.8053124
St18|St18	279.69289	-1.0473877	0.9882176	-1.0598756	0.2892012	0.7120102
Rb1cc1|Rb1cc1	12217.86115	0.3066224	0.3041177	1.0082359	0.3133412	0.7348823
4732440D04Rik|4732440D04Rik	5369.16278	-0.1463996	0.2594995	-0.5641616	0.5726441	0.9132695
Gm26901|Gm26901	163.46252	-0.7535570	1.0236447	-0.7361509	0.4616389	0.8511309
Pcmtd1|Pcmtd1	4321.05905	0.5461501	0.2827538	1.9315391	0.0534164	0.3103856
Rrs1|Rrs1	625.86245	-0.1126443	0.7061219	-0.1595253	0.8732551	0.9949371
Mybl1|Mybl1	184.22714	-2.0915101	1.2989231	-1.6101878	0.1073569	0.4493503
Adhfe1|Adhfe1	452.68051	-1.0496649	0.8432838	-1.2447350	0.2132292	0.6256462

out of 23922 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 628, 2.6% LFC < 0 (down) : 1046, 4.4% outliers [1] : 2272, 9.5% low counts [2] : 2713, 11% (mean count < 6) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

1.1.3 Craving vs Maintenance

Craving vs Maintenance CIVSA Results Overview

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
Mrpl15|Mrpl15	2162.26767	-0.2627261	0.3466449	-0.7579115	0.4485040	0.8651323
Gm39587|Gm39587	20.30129	-0.3516955	3.3976607	-0.1035111	0.9175574	0.9939251
Lypla1|Lypla1	994.90242	-0.0715356	0.8799222	-0.0812976	0.9352053	0.9951721
Xkr4|Xkr4	4459.42848	-0.6679818	0.2714805	-2.4605151	0.0138738	0.1495170
Gm39585|Gm39585	5797.38480	-1.0936780	0.2412810	-4.5327985	0.0000058	0.0002749
Tcea1|Tcea1	2657.79689	-0.5356827	0.3330620	-1.6083574	0.1077569	0.4613839
Rgs20|Rgs20	422.79650	-0.5371710	0.7618243	-0.7051114	0.4807409	0.8833096
Gm16041|Gm16041	21.30106	0.4549475	3.4906011	0.1303350	0.8963014	0.9909408
Atp6v1h|Atp6v1h	10659.10495	-0.6141458	0.2967693	-2.0694383	0.0385050	0.2733236
LOC108167788|LOC108167788	11.33585	-0.2865007	3.5020543	-0.0818093	0.9347983	0.9951721
Npbwr1|Npbwr1	924.87540	0.3369643	0.6640179	0.5074627	0.6118302	0.9372345
Oprk1|Oprk1	15661.74234	-0.4369729	0.3231413	-1.3522657	0.1762903	0.5932565
St18|St18	279.69289	-0.9608973	0.9259208	-1.0377748	0.2993749	0.7560424
Rb1cc1|Rb1cc1	12217.86115	-0.3863013	0.2844705	-1.3579665	0.1744743	0.5910084
4732440D04Rik|4732440D04Rik	5369.16278	-0.3454179	0.2427896	-1.4227046	0.1548218	0.5589467
Gm26901|Gm26901	163.46252	0.5242196	0.9577914	0.5473213	0.5841580	0.9274722
Pcmtd1|Pcmtd1	4321.05905	-0.3426769	0.2644209	-1.2959526	0.1949918	0.6240924
Rrs1|Rrs1	625.86245	-0.3724096	0.6606702	-0.5636845	0.5729689	0.9244236
Mybl1|Mybl1	184.22714	3.5366832	1.2149336	2.9110094	0.0036026	0.0603776
Adhfe1|Adhfe1	452.68051	0.8959769	0.7889411	1.1356702	0.2560946	0.7074992

out of 23922 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 1084, 4.5% LFC < 0 (down) : 324, 1.4% outliers [1] : 2271, 9.5% low counts [2] : 2713, 11% (mean count < 6) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

1.2 Annotate Results Files and Pull Significant DEGs

This subsection is for annotation the results files to include all necessary information for downstream analyses.

After annotation I subset the results to only view the significantly differentially expressed genes with an adjusted p-value > 0.05.

## Subset for significant genes only
results_sig <- subset(res, padj < 0.05)
results_sig2 <- subset(res2, padj < 0.05)
results_sig3 <- subset(res3, padj < 0.05)

1.2.1 Craving vs FoodTrained

Significant DEGs Craving vs FoodTrained

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	symbol	description	entrez	ensembl
Gm28836|Gm28836	24.57062	18.512895	3.6825760	5.027159	0.0000005	0.0000748	Gm28836	predicted gene 28836	102635593	NA
Vwa3b|Vwa3b	68.09876	9.322079	2.2726944	4.101774	0.0000410	0.0048096	Vwa3b	von Willebrand factor A domain containing 3B	70853	ENSMUSG0….
4930439A04Rik|4930439A04Rik	239.71050	6.086391	1.0543294	5.772760	0.0000000	0.0000025	4930439A04Rik	RIKEN cDNA 4930439A04 gene	78119	NA
LOC108167622|LOC108167622	55.54376	21.271643	2.7752955	7.664641	0.0000000	0.0000000	LOC108167622	NA	NA
Pard3bos2|Pard3bos2	48.12983	-22.863164	3.6323379	-6.294338	0.0000000	0.0000002	Pard3bos2	par-3 family cell polarity regulator beta, opposite strand 2	105243940	ENSMUSG0….
Col4a3|Col4a3	29.96706	19.601203	3.5629545	5.501390	0.0000000	0.0000079	Col4a3	collagen, type IV, alpha 3	12828	ENSMUSG0….
Dnajb3|Dnajb3	70.36681	9.093523	2.0877638	4.355628	0.0000133	0.0017489	Dnajb3	DnaJ heat shock protein family (Hsp40) member B3	15504	ENSMUSG0….
Asb18|Asb18	12.85015	18.822115	3.6835623	5.109759	0.0000003	0.0000497	Asb18	ankyrin repeat and SOCS box-containing 18	208372	ENSMUSG0….
Gm41920|Gm41920	91.39199	20.766926	3.6195797	5.737386	0.0000000	0.0000030	Gm41920	predicted gene, 41920	105246666	NA
LOC105246710|LOC105246710	41.29964	21.372039	3.6789651	5.809253	0.0000000	0.0000022	LOC105246710	NA	NA
LOC108167732|LOC108167732	59.78860	21.982685	3.2366754	6.791749	0.0000000	0.0000000	LOC108167732	NA	NA
Fam129a|Fam129a	58.26340	9.542099	2.6179113	3.644928	0.0002675	0.0256134	Fam129a	NA	NA
Fmo1|Fmo1	51.26689	21.615385	2.7780474	7.780783	0.0000000	0.0000000	Fmo1	flavin containing monooxygenase 1	14261	ENSMUSG0….
Gorab|Gorab	666.84329	2.774594	0.5327144	5.208407	0.0000002	0.0000307	Gorab	golgin, RAB6-interacting	98376	ENSMUSG0….
Dnah14|Dnah14	41.60299	21.334681	3.4180513	6.241767	0.0000000	0.0000003	Dnah14	dynein, axonemal, heavy chain 14	240960	ENSMUSG0….
Atf3|Atf3	433.56152	-2.808330	0.7372279	-3.809311	0.0001394	0.0145805	Atf3	activating transcription factor 3	11910	ENSMUSG0….
B230208H11Rik|B230208H11Rik	48.28040	21.761939	2.7106546	8.028297	0.0000000	0.0000000	B230208H11Rik	RIKEN cDNA B230208H11 gene	320273	ENSMUSG0….
LOC108167797|LOC108167797	82.55280	9.456510	2.0988898	4.505482	0.0000066	0.0009335	LOC108167797	NA	NA
Fam26f|Fam26f	45.72461	21.035286	2.7780345	7.572003	0.0000000	0.0000000	Fam26f	NA	NA
Mical1|Mical1	109.45019	10.585487	2.6621850	3.976240	0.0000700	0.0077657	Mical1	microtubule associated monooxygenase, calponin and LIM domain containing 1	171580	ENSMUSG0….

1.2.2 Maintenance vs FoodTrained

Significant DEGs Maintenance vs FoodTrained

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	symbol	description	entrez	ensembl
Gm39585|Gm39585	5797.38480	1.179803	0.2580070	4.572755	0.0000048	0.0001973	Gm39585	predicted gene, 39585	105243853	NA
Snhg6|Snhg6	522.81858	1.260214	0.4234933	2.975758	0.0029227	0.0446701	Snhg6	small nucleolar RNA host gene 6	73824	ENSMUSG0….
Trpa1|Trpa1	47.61966	-22.462571	3.6001396	-6.239361	0.0000000	0.0000001	Trpa1	transient receptor potential cation channel, subfamily A, member 1	277328	ENSMUSG0….
Rpl7|Rpl7	2486.96588	1.005240	0.3038678	3.308150	0.0009391	0.0184408	Rpl7	ribosomal protein L7	19989	ENSMUSG0….
Gm28836|Gm28836	24.57062	20.693225	3.9082681	5.294730	0.0000001	0.0000064	Gm28836	predicted gene 28836	102635593	NA
Gm32790|Gm32790	44.07127	9.563380	2.2173604	4.312957	0.0000161	0.0005913	Gm32790	predicted gene, 32790	102635458	NA
Fer1l5|Fer1l5	230.33092	-2.687277	0.8236933	-3.262473	0.0011044	0.0209149	Fer1l5	fer-1-like 5 (C. elegans)	100534273	ENSMUSG0….
Sema4c|Sema4c	1298.41081	-1.822125	0.6102459	-2.985887	0.0028276	0.0434624	Sema4c	sema domain, immunoglobulin domain (Ig), transmembrane domain (TM) and short cytoplasmic domain, (semaphorin) 4C	20353	ENSMUSG0….
Vwa3b|Vwa3b	68.09876	8.717777	2.3896076	3.648204	0.0002641	0.0065801	Vwa3b	von Willebrand factor A domain containing 3B	70853	ENSMUSG0….
4930439A04Rik|4930439A04Rik	239.71050	7.561383	1.1113671	6.803677	0.0000000	0.0000000	4930439A04Rik	RIKEN cDNA 4930439A04 gene	78119	NA
Il1r2|Il1r2	73.49166	9.935852	2.8549805	3.480182	0.0005011	0.0112029	Il1r2	interleukin 1 receptor, type II	16178	ENSMUSG0….
Il1rl2|Il1rl2	153.74171	-24.241136	3.7848523	-6.404777	0.0000000	0.0000000	Il1rl2	interleukin 1 receptor-like 2	107527	ENSMUSG0….
LOC108167622|LOC108167622	55.54376	21.760261	2.9333442	7.418244	0.0000000	0.0000000	LOC108167622	NA	NA
Gm10561|Gm10561	31.31746	-8.823060	2.8132289	-3.136275	0.0017111	0.0295377	Gm10561	predicted gene 10561	628004	NA
Clk1|Clk1	11902.12420	1.056339	0.3336618	3.165897	0.0015461	0.0272603	Clk1	CDC-like kinase 1	12747	ENSMUSG0….
Stradb|Stradb	1233.07664	1.338207	0.3402418	3.933107	0.0000839	0.0024506	Stradb	STE20-related kinase adaptor beta	227154	ENSMUSG0….
Fam117b|Fam117b	1077.36822	1.740967	0.5934892	2.933444	0.0033522	0.0491723	Fam117b	family with sequence similarity 117, member B	72750	ENSMUSG0….
D230017M19Rik|D230017M19Rik	353.30341	1.817985	0.5121487	3.549721	0.0003856	0.0090159	D230017M19Rik	RIKEN cDNA D230017M19 gene	320933	NA
Arpc2|Arpc2	3532.07869	-1.115469	0.2391510	-4.664286	0.0000031	0.0001315	Arpc2	actin related protein 2/3 complex, subunit 2	76709	ENSMUSG0….
Cnppd1|Cnppd1	1597.07545	-1.205071	0.3232300	-3.728214	0.0001928	0.0049685	Cnppd1	cyclin Pas1/PHO80 domain containing 1	69171	ENSMUSG0….

1.2.3 Craving vs Maintenance

Significant Craving vs Maintenance

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj	symbol	description	entrez	ensembl
Gm39585|Gm39585	5797.38480	-1.0936780	0.2412810	-4.532799	0.0000058	0.0002749	Gm39585	predicted gene, 39585	105243853	NA
Trpa1|Trpa1	47.61966	20.0711689	3.3932805	5.914975	0.0000000	0.0000004	Trpa1	transient receptor potential cation channel, subfamily A, member 1	277328	ENSMUSG0….
Rpl7|Rpl7	2486.96588	-1.0763933	0.2841797	-3.787721	0.0001520	0.0049557	Rpl7	ribosomal protein L7	19989	ENSMUSG0….
1700001G17Rik|1700001G17Rik	704.98525	1.5933623	0.4255143	3.744557	0.0001807	0.0056755	1700001G17Rik	RIKEN cDNA 1700001G17 gene	67503	ENSMUSG0….
Fer1l5|Fer1l5	230.33092	2.7365069	0.7714625	3.547168	0.0003894	0.0105650	Fer1l5	fer-1-like 5 (C. elegans)	100534273	ENSMUSG0….
Sema4c|Sema4c	1298.41081	1.9826591	0.5709143	3.472779	0.0005151	0.0131477	Sema4c	sema domain, immunoglobulin domain (Ig), transmembrane domain (TM) and short cytoplasmic domain, (semaphorin) 4C	20353	ENSMUSG0….
Il1rl2|Il1rl2	153.74171	22.2736364	3.5649386	6.247972	0.0000000	0.0000001	Il1rl2	interleukin 1 receptor-like 2	107527	ENSMUSG0….
Fam126b|Fam126b	8348.97120	-0.8264248	0.2373035	-3.482565	0.0004966	0.0127789	Fam126b	NA	NA
Pard3bos2|Pard3bos2	48.12983	-23.7390363	3.6316177	-6.536766	0.0000000	0.0000000	Pard3bos2	par-3 family cell polarity regulator beta, opposite strand 2	105243940	ENSMUSG0….
Zdbf2|Zdbf2	3291.47772	-0.9550029	0.3014901	-3.167610	0.0015370	0.0313656	Zdbf2	zinc finger, DBF-type containing 2	73884	ENSMUSG0….
Klf7|Klf7	3114.99312	-0.6650793	0.2115799	-3.143396	0.0016700	0.0335025	Klf7	Kruppel-like factor 7 (ubiquitous)	93691	ENSMUSG0….
Map2|Map2	15113.61849	-0.6774675	0.1781094	-3.803658	0.0001426	0.0047204	Map2	microtubule-associated protein 2	17756	ENSMUSG0….
Gm31047|Gm31047	43.06422	19.1697461	3.6817645	5.206674	0.0000002	0.0000130	Gm31047	predicted gene, 31047	102633145	NA
Arpc2|Arpc2	3532.07869	0.7956339	0.2237915	3.555246	0.0003776	0.0103047	Arpc2	actin related protein 2/3 complex, subunit 2	76709	ENSMUSG0….
Cnppd1|Cnppd1	1597.07545	0.9103462	0.3024970	3.009439	0.0026173	0.0473299	Cnppd1	cyclin Pas1/PHO80 domain containing 1	69171	ENSMUSG0….
Dnajb2|Dnajb2	1487.62355	1.1910673	0.3095799	3.847366	0.0001194	0.0040449	Dnajb2	DnaJ heat shock protein family (Hsp40) member B2	56812	ENSMUSG0….
Atg9a|Atg9a	8813.08524	2.0179237	0.4649911	4.339704	0.0000143	0.0006226	Atg9a	autophagy related 9A	245860	ENSMUSG0….
Resp18|Resp18	56768.10866	-0.8522294	0.2302091	-3.701980	0.0002139	0.0065238	Resp18	regulated endocrine-specific protein 18	19711	ENSMUSG0….
Chpf|Chpf	11304.97381	2.0829961	0.3845635	5.416520	0.0000001	0.0000049	Chpf	chondroitin polymerizing factor	74241	ENSMUSG0….
Tmem198|Tmem198	3311.19558	1.5407766	0.3702558	4.161384	0.0000316	0.0012693	Tmem198	transmembrane protein 198	319998	ENSMUSG0….

1.3 Visualize Shrinkage Estimations

Following subsetting, I want to see how the data looks based on effect size so I have run a shrinkage model using the “apeglm” method which normalizes counts based on effect size without major changes to data structure. I have included non-shrunk models as well for comparison.

1.3.1 Non-Shrunken L2FC

1.3.1.1 Craving vs FoodTrained

Past versions of "Non-Shrunken Estimate Plots-1.png"

Version	Author	Date
dc04b15	avm27	2024-07-17

1.3.1.2 Maintenance vs FoodTrained

Past versions of "Non-Shrunken Estimate Plots2-1.png"

Version	Author	Date
dc04b15	avm27	2024-07-17

1.3.1.3 Craving vs Maintenance

Past versions of "Non-Shrunken Estimate Plots3-1.png"

Version	Author	Date
dc04b15	avm27	2024-07-17

1.3.2 Shrunken L2FC

1.3.2.1 Craving vs FoodTrained

Past versions of "Shrunken Estimate Plots-1.png"

Version	Author	Date
dc04b15	avm27	2024-07-17

1.3.2.2 Maintenance vs FoodTrained

Past versions of "Shrunken Estimate Plots2-1.png"

Version	Author	Date
dc04b15	avm27	2024-07-17

1.3.2.3 Craving vs Maintenance

Past versions of "Shrunken Estimate Plots3-1.png"

Version	Author	Date
dc04b15	avm27	2024-07-17

2 Looking at Differentially Expressed Genes

We can generate a matrix to identify where our differentially expressed genes lie between each group based on an alpha value of 0.05.

## Matrix Visualizaattion
x <- vsDEGMatrix(
  data = dds, padj = 0.05, d.factor = "condition",
  type = "deseq", title = TRUE, legend = TRUE, grid = TRUE
)

Past versions of "Matrix of Differentially Expressed genes across groups-1.png"

Version	Author	Date
dc04b15	avm27	2024-07-17

3 Data Visualization

3.1 Normalized Counts and Plots of Genes of Interest

level_order <- c("FoodTrained", "Maintenance", "Craving") # this vector might be useful for other plots/analyses

normalized_counts <- counts(dds, normalized = TRUE)

3.1.1 Grin1

#### Grin1

d2 <- plotCounts(dds,
  gene = "Grin1|Grin1", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Grin1 Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Past versions of "Selected Genes-1.png"

Version	Author	Date
12fbe39	avm27	2024-07-17
25a76da	avm27	2024-07-17
cb49269	avm27	2024-07-17
21b70e3	avm27	2024-07-17
a1c090c	avm27	2024-07-17
dc04b15	avm27	2024-07-17

3.1.2 Gabbr1

#### Gabbr1


d2 <- plotCounts(dds,
  gene = "Gabbr1|Gabbr1", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Gabbr1 Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Past versions of "Selected Genes2-1.png"

Version	Author	Date
12fbe39	avm27	2024-07-17
25a76da	avm27	2024-07-17
cb49269	avm27	2024-07-17
21b70e3	avm27	2024-07-17
a1c090c	avm27	2024-07-17
dc04b15	avm27	2024-07-17

3.1.3 Th

#### Th


d2 <- plotCounts(dds,
  gene = "Th|Th", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Th Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Past versions of "Selected Genes3-1.png"

Version	Author	Date
12fbe39	avm27	2024-07-17
25a76da	avm27	2024-07-17
cb49269	avm27	2024-07-17
21b70e3	avm27	2024-07-17
a1c090c	avm27	2024-07-17
dc04b15	avm27	2024-07-17

3.1.4 Ehmt2

#### Ehmt2

d2 <- plotCounts(dds,
  gene = "Ehmt2|Ehmt2", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Ehmt2 Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Past versions of "Selected Genes4-1.png"

Version	Author	Date
12fbe39	avm27	2024-07-17
25a76da	avm27	2024-07-17
cb49269	avm27	2024-07-17
21b70e3	avm27	2024-07-17
a1c090c	avm27	2024-07-17
dc04b15	avm27	2024-07-17

3.1.5 Setdb1

#### Setdb1

d2 <- plotCounts(dds,
  gene = "Setdb1|Setdb1", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Setdb1 Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Past versions of "Selected Genes5-1.png"

Version	Author	Date
12fbe39	avm27	2024-07-17
25a76da	avm27	2024-07-17
cb49269	avm27	2024-07-17
21b70e3	avm27	2024-07-17
a1c090c	avm27	2024-07-17
dc04b15	avm27	2024-07-17

3.1.6 Atf7ip

#### Atf7ip

d2 <- plotCounts(dds,
  gene = "Atf7ip|Atf7ip", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Atf7ip Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Past versions of "Selected Genes6-1.png"

Version	Author	Date
12fbe39	avm27	2024-07-17
25a76da	avm27	2024-07-17
cb49269	avm27	2024-07-17
21b70e3	avm27	2024-07-17
a1c090c	avm27	2024-07-17
dc04b15	avm27	2024-07-17

3.2 Principal Components Analysis (PCA)

For visualization I am converting all analyses to logarithmic. Following this I have completed a PCA for each analyses based on all genes together.

Past versions of "Plot PCA-1.png"

Version	Author	Date
12fbe39	avm27	2024-07-17
25a76da	avm27	2024-07-17
cb49269	avm27	2024-07-17
21b70e3	avm27	2024-07-17
a1c090c	avm27	2024-07-17
dc04b15	avm27	2024-07-17

3.3 Venn Diagrams

3.3.1 Up-Regulated Genes

myCol <- brewer.pal(3, "Pastel2")

results_sig_Table_up <- subset(results_sig, log2FoldChange > 0)
results_sig2_Table_up <- subset(results_sig2, log2FoldChange > 0)
results_sig3_Table_up <- subset(results_sig3, log2FoldChange > 0)

new_ressig1_up  <- rownames(results_sig_Table_up)
new_ressig2_up <- rownames(results_sig2_Table_up)
new_ressig3_up  <- rownames(results_sig3_Table_up)

# Generate Venn Diagram 
venn.plot.up <- venn.diagram(
  x = list(
    "Crv_vs_FT" = new_ressig1_up,
    "Maint_vs_FT" = new_ressig2_up,
    "Crv_vs_Maint" = new_ressig3_up),
    fill = myCol,
    cex = .8,
    fontface = "bold",
    fontfamily = "sans",
    cat.cex = 0.8,
    cat.fontface = "bold",
    cat.default.pos = "outer",
    cat.pos = c(0, 0, 0),
    cat.fontfamily = "sans",
    rotation = 1,
  category.names = c("Craving_vs_FoodTrained", "Maintenance_vs_FoodTrained", "Craving_vs_Maintenance"),
  filename = NULL,
  output = TRUE
)

# Save Venn Diagram
#pdf("plots/CIVSA_venn_diagram_Up.pdf")
grid.draw(venn.plot.up)

Past versions of unnamed-chunk-1-1.png

Version	Author	Date
dc04b15	avm27	2024-07-17

#dev.off()

# Display Venn Diagram in RStudio
#grid.draw(venn.plot.up)

3.3.2 Down-Regulated Genes

results_sig_Table_down <- subset(results_sig, log2FoldChange < 0)
results_sig2_Table_down <- subset(results_sig2, log2FoldChange < 0)
results_sig3_Table_down <- subset(results_sig3, log2FoldChange < 0)

new_ressig1_down  <- rownames(results_sig_Table_down)
new_ressig2_down  <- rownames(results_sig2_Table_down)
new_ressig3_down  <- rownames(results_sig3_Table_down)


# Generate Venn Diagram UP
venn.plot.down <- venn.diagram(
  x = list(
    "Crv_vs_FT" = new_ressig1_down,
    "Maint_vs_FT" = new_ressig2_down,
    "Crv_vs_Maint" = new_ressig3_down),
    fill = myCol,
    cex = .8,
    fontface = "bold",
    fontfamily = "sans",
    cat.cex = 0.8,
    cat.fontface = "bold",
    cat.default.pos = "outer",
    cat.pos = c(0, 0, 0),
    cat.fontfamily = "sans",
    rotation = 1,
  category.names = c("Craving_vs_FoodTrained", "Maintenance_vs_FoodTrained", "Craving_vs_Maintenance"),
  filename = NULL,
  output = TRUE
)

# Save Venn Diagram
#pdf("plots/CIVSA_venn_diagram_Down.pdf")
#grid.draw(venn.plot.down)
#dev.off()

# Display Venn Diagram in RStudio
grid.draw(venn.plot.down)

Past versions of unnamed-chunk-2-1.png

Version	Author	Date
dc04b15	avm27	2024-07-17

3.4 Heatmaps of Differentially Expressed Genes

ressig1 <- as.data.frame(results_sig)
ressig2 <- as.data.frame(results_sig2)
ressig3 <- as.data.frame(results_sig3)

ressig1$Comparison <- "Craving vs FoodTrained"
ressig2$Comparison <- "Maintenance vs FoodTrained"
ressig3$Comparison <- "Craving vs Maintenance"

ressig1_highlfc_pos <- ressig1[ressig1$log2FoldChange >= 1.5,]
ressig1_highlfc_neg <- ressig1[ressig1$log2FoldChange <= -1.5,]
ressig2_highlfc_pos <- ressig2[ressig2$log2FoldChange >= 1.5,]
ressig2_highlfc_neg <- ressig2[ressig2$log2FoldChange <= -1.5,]
ressig3_highlfc_pos <- ressig3[ressig3$log2FoldChange >= 1.5,]
ressig3_highlfc_neg <- ressig3[ressig3$log2FoldChange <= -1.5,]

temp_sigdegshm <- rbind(ressig1_highlfc_pos, ressig1_highlfc_neg)
temp_sigdegshm2 <- rbind(temp_sigdegshm, ressig2_highlfc_pos)
temp_sigdegshm3 <- rbind(temp_sigdegshm2, ressig2_highlfc_neg)
temp_sigdegshm4 <- rbind(temp_sigdegshm3, ressig3_highlfc_pos)
temp_sigdegshm5 <- rbind(temp_sigdegshm4, ressig3_highlfc_neg)

temp_sigdegshm6 <- temp_sigdegshm5[temp_sigdegshm5$lfcSE <=1,]

dup_rows <- duplicated(temp_sigdegshm6$symbol)
data_dedup <- temp_sigdegshm6[!dup_rows, ]

rows <- rownames(data_dedup)

select <- order(rowMeans(counts(dds,normalized=TRUE)),
                decreasing=TRUE)#[1:1500]

mat <- assay(ddsMat_rlog)[select,]

data_mod <- mat[rownames(mat) %in% rows, ] 

matrownames_data_mod <- data.frame(rownames(data_mod))

matrownames_data_mod <- matrownames_data_mod %>% separate(rownames.data_mod., c("id", "symbol"),sep ="([|])")

samples <- as.tibble(read.csv("data/sample_info_CIVSA.csv", sep = ","))

annot_col <- samples %>%
  column_to_rownames('samplename') %>%
  as.data.frame()

ann_colors = list(
  condition = c(FoodTrained = "#0000FF", Maintenance = "#AD07E3", Craving = "#FF8000"))

3.4.1 Unsupervised Clustering

pheatmap::pheatmap(
  mat = data_mod,
  color = colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(155),
  scale = "row",
  annotation_colors = ann_colors,
  annotation_col = annot_col,
  annotation_legend = T,
  fontsize = 10, 
  show_colnames = T,
  cluster_rows = T,
  cluster_cols = T,
  show_rownames = F,
  clustering_distance_rows = "correlation",
  clustering_method = "ward.D",
  main = "Heatmap of Significant DEGs -- Unsupervised Clustering",
  cluster_row_slices = TRUE,
  show_row_dend = FALSE,
  cutree_rows = 4,
  cutree_cols = 3,
  treeheight_row = 10,
  treeheight_col = 10,
  #filename = "plots/CIVSA_Heatmap_DEGs_Final_Condensed_Unsupervised.pdf",
  legend_labels = c("FoodTrained", "Maintenance", "Craving"),
)

Past versions of "Filter Heatmaps by DEGs2 Unsupervised-1.png"

Version	Author	Date
12fbe39	avm27	2024-07-17
25a76da	avm27	2024-07-17
21b70e3	avm27	2024-07-17
a1c090c	avm27	2024-07-17
dc04b15	avm27	2024-07-17

3.4.2 Supervised Clustering

pheatmap::pheatmap(
  mat = data_mod,
  color = colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(155),
  scale = "row",
  annotation_col = annot_col,
  annotation_colors = ann_colors,
  fontsize = 10,
  show_colnames = T,
  cluster_rows = T,
  cluster_cols = F,
  annotation_legend = T,
  legend = T,
  show_rownames = F,
  clustering_distance_rows = "correlation",
  clustering_method = "ward.D",
  main = "Heatmap of Significant DEGs -- Supervised Clustering",
  cluster_row_slices = TRUE,
  show_row_dend = FALSE,
  cutree_rows = 4,
  cutree_cols = 3,
  gaps_col = c(3, 6),
  treeheight_row = 10,
  treeheight_col = 10,
  #filename = "plots/CIVSA_Heatmap_DEGs_Final_Condensed_Supervised.pdf",
  legend_labels = c("FoodTrained", "Maintenance", "Craving"),
)

Past versions of "Filter Heatmaps by DEGs2 Condensed-Supervised-1.png"

Version	Author	Date
12fbe39	avm27	2024-07-17
25a76da	avm27	2024-07-17
21b70e3	avm27	2024-07-17
a1c090c	avm27	2024-07-17
dc04b15	avm27	2024-07-17

3.5 Volcano Plots

3.5.1 Visualize Results: Volcano Plots (Non-Interactive)

Volcano plots are generated by gathering the FDR corrected p-values from each analysis. The adjusted p-values undergo a -log10 transformation to generate these FDR corrected values. Any rows containing “NA” are ommitted from analysis. Data points are colored based on increasing or decreasing value and plotted using ggplot2. Adjusted p-value cutoff is 1.3 following log transformation.(-log10 0.05 = 1.3)

Past versions of "Volcano Plots-1.png"

Version	Author	Date
a1c090c	avm27	2024-07-17

3.5.2 Visualize Results: Volcano Plots (Interactive)

Interactive Volcano plots are generated by Using the same method as above, but are plotted so that genes may be identified and log2FC can be seen from each gene. Plots can be used to customize which genes appear based on what you would like to see.

4 Pathway Analysis

4.1 Pathway Analysis GO

First, the background for the GO universe must be set. In order to do so we will subset the first results file “res” for any genes that map to entrez IDs and then we only take the unique ones. This leaves us with ~31,000 genes for the remaining background which we list as our midbrain dopamine neuronal expressed genes.

Pathway analysis is conducted using enrichGo and enrichKEGG based on a p-value of 0.05 and a q-value of 0.10 and including all significantly differentially expressed genes.

# Remove any genes that have greater than 1.5 LFC Standard Error
results_sig_entrez <- subset(results_sig, lfcSE <= 1.5)
results_sig_entrez2 <- subset(results_sig2, lfcSE <= 1.5)
results_sig_entrez3 <- subset(results_sig3, lfcSE <= 1.5)

# Remove any genes that have greater than adj.pvalue 0.05
results_sig_entrez <- subset(results_sig_entrez, padj <= 0.05)
results_sig_entrez2 <- subset(results_sig_entrez2, padj <= 0.05)
results_sig_entrez3 <- subset(results_sig_entrez3, padj <= 0.05)

# Remove any genes that do not have any entrez identifiers
results_sig_entrez <- subset(results_sig_entrez, is.na(entrez) == FALSE)
results_sig_entrez2 <- subset(results_sig_entrez2, is.na(entrez) == FALSE)
results_sig_entrez3 <- subset(results_sig_entrez3, is.na(entrez) == FALSE)

# Create a matrix of gene log2 fold changes
gene_matrix <- results_sig_entrez$log2FoldChange
gene_matrix2 <- results_sig_entrez2$log2FoldChange
gene_matrix3 <- results_sig_entrez3$log2FoldChange

# Add the entrezID's as names for each logFC entry
names(gene_matrix) <- results_sig_entrez$entrez
names(gene_matrix2) <- results_sig_entrez2$entrez
names(gene_matrix3) <- results_sig_entrez3$entrez


go <- enrichGO(
  gene = names(gene_matrix),
  OrgDb = "org.Mm.eg.db",
  readable = T,
  universe = gene_matrix_background,
  ont = "ALL",
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.10
)

go2 <- enrichGO(
  gene = names(gene_matrix2),
  OrgDb = "org.Mm.eg.db",
  readable = T,
  universe = gene_matrix_background,
  ont = "ALL",
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.10
)

go3 <- enrichGO(
  gene = names(gene_matrix3),
  OrgDb = "org.Mm.eg.db",
  readable = T,
  universe = gene_matrix_background,
  ont = "ALL",
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.10
)

4.1.1 Craving vs FoodTrained

kable(head(go, 10), caption = "Top 10 GO Pathways for Craving vs FoodTrained", format = "html", align = "l", booktabs = TRUE) %>%
  kable_styling(html_font = "Arial", font_size = 10, bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
  scroll_box(width = "1000px", height = "500px")

Top 10 GO Pathways for Craving vs FoodTrained

	ONTOLOGY	ID	Description	GeneRatio	BgRatio	pvalue	p.adjust	qvalue	geneID	Count
GO:0005344	MF	GO:0005344	oxygen carrier activity	2/50	11/20962	0.0003025	0.0208720	0.0168759	Hba-a1/Hbb-bt	2
GO:0046974	MF	GO:0046974	histone methyltransferase activity (H3-K9 specific)	2/50	11/20962	0.0003025	0.0208720	0.0168759	Ash1l/Setdb1	2
GO:0016209	MF	GO:0016209	antioxidant activity	3/50	75/20962	0.0007641	0.0329752	0.0266618	Hba-a1/Nxn/Hbb-bt	3
GO:0008170	MF	GO:0008170	N-methyltransferase activity	3/50	81/20962	0.0009558	0.0329752	0.0266618	Pemt/Ash1l/Setdb1	3
GO:0019825	MF	GO:0019825	oxygen binding	2/50	23/20962	0.0013662	0.0377084	0.0304889	Hba-a1/Hbb-bt	2
GO:0051087	MF	GO:0051087	chaperone binding	3/50	106/20962	0.0020727	0.0476730	0.0385457	Sdf2l1/Tfrc/H2-K1	3

4.1.2 Maintenance vs FoodTrained

kable(head(go2, 10), caption = "Top 10 GO Pathways for Maintenance vs FoodTrained", format = "html", align = "l", booktabs = TRUE) %>%
  kable_styling(html_font = "Arial", font_size = 10, bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
  scroll_box(width = "1000px", height = "500px")

Top 10 GO Pathways for Maintenance vs FoodTrained

	ONTOLOGY	ID	Description	GeneRatio	BgRatio	pvalue	p.adjust	qvalue	geneID	Count
GO:0031346	BP	GO:0031346	positive regulation of cell projection organization	36/660	429/21359	1.00e-07	0.0002952	0.0002658	Arpc2/Obsl1/Cd24a/Palm/Cobl/Dlg4/Itga3/Mapt/Sphk1/Hnrnpk/Plxnb2/Twf1/Opa1/Ptk7/Ppp2r5b/Caprin1/Src/Mfn1/Ptprd/Wrap73/Ppp1r35/Fscn1/Rac1/Met/Megf8/Ap2a1/Ptpn5/Rgma/Apbb1/Cx3cl1/Katnb1/Ndrg4/Trim67/Sema7a/Map2k1/Golga4	36
GO:0010976	BP	GO:0010976	positive regulation of neuron projection development	22/660	222/21359	1.60e-06	0.0036604	0.0032953	Cd24a/Cobl/Itga3/Mapt/Sphk1/Hnrnpk/Plxnb2/Twf1/Opa1/Ptk7/Ppp2r5b/Caprin1/Mfn1/Met/Ap2a1/Ptpn5/Rgma/Apbb1/Cx3cl1/Katnb1/Ndrg4/Trim67	22
GO:0098742	BP	GO:0098742	cell-cell adhesion via plasma-membrane adhesion molecules	21/660	213/21359	3.10e-06	0.0045449	0.0040916	Obsl1/Cadm3/Cd24a/Cdh24/Pcdh8/Plxnb2/Igsf11/Ptprs/Pcdha4/Mbp/Cdh22/Fat4/Ptprd/Clstn3/Nectin2/Cadm4/Lrrc4b/Ric8a/Cx3cl1/Cdh13/Pcdh19	21
GO:0046777	BP	GO:0046777	protein autophosphorylation	22/660	241/21359	6.30e-06	0.0056489	0.0050855	Clk1/Mapkapk2/Cdk12/Aatk/Calm1/Map3k1/Pim3/Acvr1b/Dyrk1a/Ddr1/Ppp2r5b/Map3k11/Stk39/Csnk2a1/Src/Jun/Epha8/Pink1/Atp13a2/Met/Calm3/Ulk3	22
GO:0071805	BP	GO:0071805	potassium ion transmembrane transport	20/660	206/21359	6.60e-06	0.0056489	0.0050855	Calm1/Kcnk9/Kcnh8/Kcnk12/Slc12a2/Stk39/Slc12a5/Kcnq2/Hcn3/Kcna2/Slc24a2/Kcnab2/Kcnh2/Abcb8/Kcnip4/Pkd2/Calm3/Slc9a5/Nedd4/Kcnd1	20
GO:0050808	BP	GO:0050808	synapse organization	34/660	490/21359	1.08e-05	0.0056489	0.0050855	Obsl1/Ngef/Palm/Dlg4/Cdk5r1/Itga3/Mapt/Rock2/Hnrnpk/Pcdh8/Ywhaz/Adgrb1/Plxnb2/Etv5/Opa1/Ptprs/Pcdhgc5/Sptbn2/Nrxn2/Caprin1/Mfn1/Bcan/Vcp/Ptprd/Ssh1/Setd5/Iqsec3/Clstn3/Lrrc4b/Ppfia3/Apbb1/Nedd4/Tpbg/Nlgn3	34
GO:0007156	BP	GO:0007156	homophilic cell adhesion via plasma membrane adhesion molecules	14/660	114/21359	1.14e-05	0.0056489	0.0050855	Obsl1/Cadm3/Cdh24/Pcdh8/Plxnb2/Igsf11/Pcdha4/Cdh22/Fat4/Clstn3/Nectin2/Cadm4/Cdh13/Pcdh19	14
GO:0018107	BP	GO:0018107	peptidyl-threonine phosphorylation	14/660	114/21359	1.14e-05	0.0056489	0.0050855	Clk1/Cdc42bpa/Dyrk2/Cdk5r1/Sphk1/Rock2/Dmtn/Acvr1b/Dyrk1a/Stk39/Csnk2a1/Dgkq/Met/Cadm4	14
GO:0032271	BP	GO:0032271	regulation of protein polymerization	19/660	196/21359	1.14e-05	0.0056489	0.0050855	Arpc2/Camsap2/Inpp5j/Mapt/Hsp90aa1/Map3k1/Dmtn/Twf1/Dyrk1a/Snx9/Prkce/Ssh1/Arpc3/Rac1/Met/Mtpn/Eml2/Rps3/Cotl1	19
GO:0043254	BP	GO:0043254	regulation of protein-containing complex assembly	30/660	410/21359	1.30e-05	0.0057749	0.0051989	Arpc2/Camsap2/Cd24a/Ncln/Inpp5j/Mapt/Hsp90aa1/Map3k1/Dmtn/Twf1/Wnt10b/Dyrk1a/Snx9/Prkce/Ppp2r5b/Src/Arhgef2/Vcp/Pink1/Ssh1/Arpc3/Stx1a/Fscn1/Rac1/Met/Mtpn/Eml2/Rps3/Prrt2/Cotl1	30

4.1.3 Craving vs Maintenance

kable(head(go3, 10), caption = "Top 10 GO Pathways for Craving vs Maintenance", format = "html", align = "l", booktabs = TRUE) %>%
  kable_styling(html_font = "Arial", font_size = 10, bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
  scroll_box(width = "1000px", height = "500px")

Top 10 GO Pathways for Craving vs Maintenance

	ONTOLOGY	ID	Description	GeneRatio	BgRatio	pvalue	p.adjust	qvalue	geneID	Count
GO:0046777	BP	GO:0046777	protein autophosphorylation	22/607	241/21359	1.70e-06	0.0041919	0.0037612	Epha4/Mapkapk2/Mink1/Cdk12/Aatk/Calm1/Pim3/Acvr1b/Map3k12/Ddr1/Calm2/Ppp2r5b/Stk39/Src/Jun/Epha8/Pink1/Atp13a2/Met/Calm3/Hspa8/Ulk3	22
GO:0071805	BP	GO:0071805	potassium ion transmembrane transport	20/607	206/21359	1.90e-06	0.0041919	0.0037612	Ywhae/Calm1/Kcnk9/Dlg1/Kcnk12/Calm2/Bin1/Stk39/Slc12a5/Hcn3/Kcna2/Kcnab2/Akap9/Kcnh2/Abcb8/Kcnip4/Calm3/Slc9a5/Nedd4/Kcnd1	20
GO:0035304	BP	GO:0035304	regulation of protein dephosphorylation	13/607	96/21359	3.20e-06	0.0047657	0.0042761	Ywhae/Ppp1r9b/Ppp1r1b/Calm1/Ppp1r16a/Hsp90ab1/Calm2/Nsmf/Ptpa/Cry2/Pink1/Styxl1/Calm3	13
GO:0007041	BP	GO:0007041	lysosomal transport	14/607	117/21359	6.00e-06	0.0057455	0.0051553	Plekhm1/Chmp7/Vps52/Bin1/Epg5/Vps51/Vps18/Snx16/Pink1/Atp13a2/Gak/Cacng7/Hspa8/Nedd4	14
GO:0016570	BP	GO:0016570	histone modification	31/607	456/21359	7.70e-06	0.0057455	0.0051553	Elk4/Usp21/Ube2b/Srebf1/Asxl2/Hdac9/Bap1/Nipbl/Crebbp/Glyr1/Phf1/Ring1/Ubr2/Nfya/Dpy30/Cxxc1/Rps6ka4/Naa40/Ddb1/Phf20/Ash1l/Pygo2/Setdb1/Pink1/Ep400/Prmt8/Apbb1/Sin3b/Kmt2a/Setd2/Hdac8	31
GO:0035303	BP	GO:0035303	regulation of dephosphorylation	15/607	136/21359	7.80e-06	0.0057455	0.0051553	Ywhae/Ppp1r9b/Ppp1r1b/Calm1/Ppp1r16a/Hsp90ab1/Calm2/Nsmf/Ptpa/Cry2/Src/Pink1/Styxl1/Ppp1r35/Calm3	15
GO:0032271	BP	GO:0032271	regulation of protein polymerization	18/607	196/21359	1.32e-05	0.0083839	0.0075226	Map2/Arpc2/Inpp5j/Hsp90aa1/Dmtn/Twf1/Dlg1/Ttbk1/Bin1/Fchsd1/Stmn2/Gba2/Akap9/Met/Mtpn/Eml2/Rps3/Pik3r2	18
GO:0031346	BP	GO:0031346	positive regulation of cell projection organization	29/607	429/21359	1.68e-05	0.0093028	0.0083471	Arpc2/Epha4/Palm/Plekhm1/Rala/Twf1/Rrn3/Ppp2r5b/Grin1/Src/Stmn2/Wrap73/Fzd1/Styxl1/Ppp1r35/Fscn1/Met/Megf8/Ap2a1/Ptpn5/Rgma/Apbb1/Mob2/Crtc1/Cx3cl1/Katnb1/Ndrg4/Trim67/Nckipsd	29
GO:0006813	BP	GO:0006813	potassium ion transport	20/607	243/21359	2.25e-05	0.0110998	0.0099595	Ywhae/Calm1/Kcnk9/Dlg1/Kcnk12/Calm2/Bin1/Stk39/Slc12a5/Hcn3/Kcna2/Kcnab2/Akap9/Kcnh2/Abcb8/Kcnip4/Calm3/Slc9a5/Nedd4/Kcnd1	20
GO:0043666	BP	GO:0043666	regulation of phosphoprotein phosphatase activity	9/607	58/21359	3.50e-05	0.0141199	0.0126693	Ppp1r9b/Ppp1r1b/Calm1/Hsp90ab1/Calm2/Ptpa/Cry2/Styxl1/Calm3	9

5 GO Pathways Visualization

5.1 GO DotPlot

Session information

sessionInfo()

R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] plotly_4.10.2                            
 [2] gridExtra_2.3                            
 [3] ggrepel_0.9.3                            
 [4] BiocParallel_1.32.6                      
 [5] conflicted_1.2.0                         
 [6] VennDiagram_1.7.3                        
 [7] futile.logger_1.4.3                      
 [8] regionReport_1.32.0                      
 [9] kableExtra_1.3.4                         
[10] vidger_1.18.0                            
[11] tables_0.9.17                            
[12] knitr_1.43                               
[13] Glimma_2.8.0                             
[14] Mus.musculus_1.3.1                       
[15] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0
[16] OrganismDbi_1.40.0                       
[17] GenomicFeatures_1.50.4                   
[18] gplots_3.1.3                             
[19] lubridate_1.9.2                          
[20] forcats_1.0.0                            
[21] stringr_1.5.1                            
[22] purrr_1.0.2                              
[23] readr_2.1.4                              
[24] tibble_3.2.1                             
[25] tidyverse_2.0.0                          
[26] statmod_1.5.0                            
[27] tweeDEseqCountData_1.36.0                
[28] edgeR_3.40.2                             
[29] limma_3.54.2                             
[30] GOSemSim_2.24.0                          
[31] apeglm_1.20.0                            
[32] ggpubr_0.6.0                             
[33] SPIA_2.50.0                              
[34] KEGGgraph_1.58.3                         
[35] ggnewscale_0.4.9                         
[36] enrichplot_1.18.4                        
[37] tidyr_1.3.1                              
[38] WGCNA_1.72-1                             
[39] fastcluster_1.2.3                        
[40] dynamicTreeCut_1.63-1                    
[41] pathview_1.38.0                          
[42] gage_2.48.0                              
[43] dplyr_1.1.4                              
[44] topGO_2.50.0                             
[45] SparseM_1.81                             
[46] graph_1.76.0                             
[47] GO.db_3.16.0                             
[48] RColorBrewer_1.1-3                       
[49] genefilter_1.80.3                        
[50] pheatmap_1.0.12                          
[51] ggsci_3.0.0                              
[52] tximport_1.26.1                          
[53] org.Mm.eg.db_3.16.0                      
[54] AnnotationDbi_1.60.2                     
[55] DOSE_3.24.2                              
[56] ReactomePA_1.42.0                        
[57] biomaRt_2.54.1                           
[58] clusterProfiler_4.6.2                    
[59] ggplot2_3.5.0.9000                       
[60] DESeq2_1.38.3                            
[61] SummarizedExperiment_1.28.0              
[62] Biobase_2.58.0                           
[63] MatrixGenerics_1.10.0                    
[64] matrixStats_1.0.0                        
[65] GenomicRanges_1.50.2                     
[66] GenomeInfoDb_1.34.9                      
[67] IRanges_2.32.0                           
[68] S4Vectors_0.36.2                         
[69] BiocGenerics_0.44.0                      
[70] gprofiler2_0.2.2                         
[71] workflowr_1.7.0                          

loaded via a namespace (and not attached):
  [1] Hmisc_5.1-0              svglite_2.1.1            ps_1.7.5                
  [4] Rsamtools_2.14.0         foreach_1.5.2            rprojroot_2.0.3         
  [7] crayon_1.5.2             MASS_7.3-60              nlme_3.1-162            
 [10] backports_1.4.1          impute_1.72.3            rlang_1.1.3             
 [13] XVector_0.38.0           HDO.db_0.99.1            callr_3.7.3             
 [16] filelock_1.0.2           rjson_0.2.21             bit64_4.0.5             
 [19] glue_1.7.0               rngtools_1.5.2           parallel_4.2.1          
 [22] processx_3.8.2           DEFormats_1.26.0         tidyselect_1.2.1        
 [25] XML_3.99-0.14            GenomicAlignments_1.34.1 xtable_1.8-4            
 [28] magrittr_2.0.3           evaluate_0.21            bibtex_0.5.1            
 [31] cli_3.6.2                zlibbioc_1.44.0          doRNG_1.8.6             
 [34] rstudioapi_0.14          whisker_0.4.1            bslib_0.5.0             
 [37] rpart_4.1.19             derfinderHelper_1.32.0   fastmatch_1.1-3         
 [40] BiocStyle_2.26.0         lambda.r_1.2.4           treeio_1.22.0           
 [43] xfun_0.39                gson_0.1.0               cluster_2.1.4           
 [46] caTools_1.18.2           tidygraph_1.2.3          KEGGREST_1.38.0         
 [49] ape_5.7-1                Biostrings_2.66.0        png_0.1-8               
 [52] reshape_0.8.9            withr_3.0.0              bitops_1.0-7            
 [55] ggforce_0.4.1            RBGL_1.74.0              plyr_1.8.8              
 [58] coda_0.19-4              bumphunter_1.40.0        pillar_1.9.0            
 [61] cachem_1.0.8             fs_1.6.2                 graphite_1.44.0         
 [64] ellipsis_0.3.2           vctrs_0.6.5              generics_0.1.3          
 [67] tools_4.2.1              foreign_0.8-84           munsell_0.5.1           
 [70] tweenr_2.0.2             fgsea_1.24.0             DelayedArray_0.24.0     
 [73] fastmap_1.1.1            compiler_4.2.1           abind_1.4-5             
 [76] httpuv_1.6.11            rtracklayer_1.58.0       GenomeInfoDbData_1.2.9  
 [79] lattice_0.21-8           utf8_1.2.4               later_1.3.1             
 [82] BiocFileCache_2.6.1      jsonlite_1.8.7           GGally_2.1.2            
 [85] scales_1.3.0             tidytree_0.4.2           carData_3.0-5           
 [88] lazyeval_0.2.2           promises_1.2.0.1         car_3.1-2               
 [91] doParallel_1.0.17        checkmate_2.2.0          rmarkdown_2.23          
 [94] cowplot_1.1.1            webshot_0.5.5            downloader_0.4          
 [97] BSgenome_1.66.3          igraph_1.5.0             survival_3.5-5          
[100] numDeriv_2016.8-1.1      yaml_2.3.7               systemfonts_1.0.4       
[103] htmltools_0.5.5          memoise_2.0.1            VariantAnnotation_1.44.1
[106] BiocIO_1.8.0             locfit_1.5-9.8           graphlayouts_1.0.0      
[109] viridisLite_0.4.2        digest_0.6.32            rappdirs_0.3.3          
[112] knitrBootstrap_1.0.2     futile.options_1.0.1     emdbook_1.3.13          
[115] RSQLite_2.3.1            yulab.utils_0.0.6        derfinder_1.32.0        
[118] data.table_1.14.8        blob_1.2.4               preprocessCore_1.60.2   
[121] labeling_0.4.3           splines_4.2.1            Formula_1.2-5           
[124] RCurl_1.98-1.12          broom_1.0.5              hms_1.1.3               
[127] colorspace_2.1-0         base64enc_0.1-3          BiocManager_1.30.21     
[130] aplot_0.1.10             nnet_7.3-19              sass_0.4.6              
[133] Rcpp_1.0.10              mvtnorm_1.2-2            fansi_1.0.6             
[136] tzdb_0.4.0               R6_2.5.1                 lifecycle_1.0.4         
[139] formatR_1.14             curl_5.0.1               ggsignif_0.6.4          
[142] jquerylib_0.1.4          Matrix_1.5-4.1           qvalue_2.30.0           
[145] org.Hs.eg.db_3.16.0      iterators_1.0.14         RefManageR_1.4.0        
[148] htmlwidgets_1.6.2        markdown_1.7             polyclip_1.10-4         
[151] crosstalk_1.2.0          shadowtext_0.1.2         timechange_0.2.0        
[154] gridGraphics_0.5-1       reactome.db_1.82.0       rvest_1.0.3             
[157] htmlTable_2.4.1          patchwork_1.2.0.9000     bdsmatrix_1.3-6         
[160] codetools_0.2-19         gtools_3.9.4             getPass_0.2-2           
[163] prettyunits_1.1.1        dbplyr_2.3.2             gtable_0.3.5            
[166] DBI_1.1.3                git2r_0.32.0             ggfun_0.1.1             
[169] httr_1.4.6               highr_0.10               KernSmooth_2.23-21      
[172] stringi_1.8.4            progress_1.2.2           reshape2_1.4.4          
[175] farver_2.1.2             annotate_1.76.0          viridis_0.6.3           
[178] Rgraphviz_2.42.0         ggtree_3.6.2             xml2_1.3.4              
[181] bbmle_1.0.25             restfulr_0.0.15          geneplotter_1.76.0      
[184] ggplotify_0.1.1          bit_4.0.5                scatterpie_0.2.1        
[187] ggraph_2.1.0             pkgconfig_2.0.3          rstatix_0.7.2           
[190] GenomicFiles_1.34.0